rm(list=ls())

#######################################################################################################################
#######################################################################################################################

library(plyr)
source("/Users/sharpe/Desktop/Young Surgeons/code/baltab_ys.R")


#########################################
#### Read in the matched ortho data. ####


mod_match_ortho <- read.csv("/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_0_ortho_matched_mod_pats.csv")
mod_match_gensurg <- read.csv("/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_0_gensurg_matched_mod_pats.csv")

trad_mod_match <- rbind.fill(mod_match_ortho, mod_match_gensurg)

#### Give patients NA for risk scores that they were not included in. ####

trad_mod_match$prop_log_pat_mod_post_2010 <- ifelse(trad_mod_match$era_mod==1 & trad_mod_match$oct_2010_plus==1, trad_mod_match$prop_log_pat_mod_post_2010, NA)
trad_mod_match$prop_prob_pat_mod_post_2010 <- ifelse(trad_mod_match$era_mod==1 & trad_mod_match$oct_2010_plus==1, trad_mod_match$prop_prob_pat_mod_post_2010, NA)
trad_mod_match$prop_log_pat_mod_pre_2010 <- ifelse(trad_mod_match$era_mod==1 & trad_mod_match$oct_2010_plus==0, trad_mod_match$prop_log_pat_mod_pre_2010, NA)
trad_mod_match$prop_prob_pat_mod_pre_2010 <- ifelse(trad_mod_match$era_mod==1 & trad_mod_match$oct_2010_plus==0, trad_mod_match$prop_prob_pat_mod_pre_2010, NA)

#### Create a combined modern risk score. ####

trad_mod_match$prop_log_pat_mod <- ifelse(trad_mod_match$era_mod==1 & trad_mod_match$oct_2010_plus==0, trad_mod_match$prop_log_pat_mod_pre_2010, trad_mod_match$prop_log_pat_mod_post_2010)
trad_mod_match$prop_prob_pat_mod <- ifelse(trad_mod_match$era_mod==1 & trad_mod_match$oct_2010_plus==0, trad_mod_match$prop_prob_pat_mod_pre_2010, trad_mod_match$prop_prob_pat_mod_post_2010)

##############################################
#### Read in the before matched datasets. #### 



new_ortho <- read.csv(file="/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_0_ortho_new_pats.csv")
exp_ortho <- read.csv(file="/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_0_ortho_exp_pats.csv")
new_gensurg <- read.csv(file="/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_0_gensurg_new_pats.csv")
exp_gensurg <- read.csv(file="/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_0_gensurg_exp_pats.csv")

new_exp <- rbind.fill(new_ortho,new_gensurg, exp_ortho, exp_gensurg)
new_exp <- subset(new_exp, era_mod==1)


new_exp$prop_log_pat_mod_post_2010 <- ifelse(new_exp$era_mod==1 & new_exp$oct_2010_plus==1, new_exp$prop_log_pat_mod_post_2010,NA)
new_exp$prop_prob_pat_mod_post_2010 <- ifelse(new_exp$era_mod==1 & new_exp$oct_2010_plus==1, new_exp$prop_prob_pat_mod_post_2010, NA)
new_exp$prop_log_pat_mod_pre_2010 <- ifelse(new_exp$era_mod==1 & new_exp$oct_2010_plus==0, new_exp$prop_log_pat_mod_pre_2010, NA)
new_exp$prop_prob_pat_mod_pre_2010 <- ifelse(new_exp$era_mod==1 & new_exp$oct_2010_plus==0, new_exp$prop_prob_pat_mod_pre_2010, NA)

new_exp$prop_log_pat_mod <- ifelse(new_exp$era_mod==1 & new_exp$oct_2010_plus==0, new_exp$prop_log_pat_mod_pre_2010, new_exp$prop_log_pat_mod_post_2010)
new_exp$prop_prob_pat_mod <- ifelse(new_exp$era_mod==1 & new_exp$oct_2010_plus==0, new_exp$prop_prob_pat_mod_pre_2010, new_exp$prop_prob_pat_mod_post_2010)

mod_unmatch <- subset(new_exp, era_mod==1)

##################################################
#### Read in the certified_by_era variable. ######
#### Merge to matched and unmatched datasets. ####

trad_gensurg_cert <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/transfer_trad_gensurg_cert.csv")
mod_gensurg_cert <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/transfer_mod_gensurg_cert.csv")
trad_ortho_cert <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/transfer_trad_ortho_cert.csv")
mod_ortho_cert <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/transfer_mod_ortho_cert.csv")

all_cert <- rbind(trad_gensurg_cert, mod_gensurg_cert, trad_ortho_cert, mod_ortho_cert)

new_exp <- merge(x=new_exp, y=all_cert, by=c("newID_1","newID_10"), all.x=T)
trad_mod_match <- merge(x=trad_mod_match, y=all_cert, by=c("newID_1","newID_10"), all.x=T)

###########################################################################################################
#### Read in the risk and comorbidity risk scores from the level 3 match and apply them to this match. ####


risk_ortho_post_2010 <- read.csv("/Users/sharpe/Desktop/Young Surgeons/table/risk_model_mod_ortho_post_2010_11.8.2018_exact_procs.csv")
risk_ortho_pre_2010 <- read.csv("/Users/sharpe/Desktop/Young Surgeons/table/risk_model_mod_ortho_post_2010_11.8.2018_exact_procs.csv")
risk_gensurg_post_2010 <- read.csv("/Users/sharpe/Desktop/Young Surgeons/table/risk_model_mod_gensurg_post_2010_11.2.2018.csv")
risk_gensurg_pre_2010 <- read.csv("/Users/sharpe/Desktop/Young Surgeons/table/risk_model_mod_gensurg_pre_2010_11.2.2018.csv")

crisk_ortho_post_2010 <- read.csv("/Users/sharpe/Desktop/Young Surgeons/table/crisk_model_mod_ortho_post_2010_11.8.2018_exact_procs.csv")
crisk_ortho_pre_2010 <- read.csv("/Users/sharpe/Desktop/Young Surgeons/table/crisk_model_mod_ortho_post_2010_11.8.2018_exact_procs.csv")
crisk_gensurg_post_2010 <- read.csv("/Users/sharpe/Desktop/Young Surgeons/table/crisk_model_mod_gensurg_post_2010_11.2.2018.csv")
crisk_gensurg_pre_2010 <- read.csv("/Users/sharpe/Desktop/Young Surgeons/table/crisk_model_mod_gensurg_pre_2010_11.2.2018.csv")

apply_model <- function(model,dataset,varname1,varname2,era){
  
  
  #### Subset data to only the variables that are in the model. ####
  get_model_vars <- as.character(model[(2:nrow(model)),1])
  temp_data <- as.matrix(subset(dataset, select=get_model_vars))
  
  #### Save the model intercept and coefficients. ####
  intercept <- model[1,2]
  coefs <- as.matrix(model[(2:nrow(model)),2])
  
  #### Multiply the coefficients by the data. ####
  scores <- as.data.frame(temp_data %*% coefs + intercept)
  dataset <- cbind(as.data.frame(dataset), scores)
  colnames(dataset)[dim(dataset)[2]] <- varname1
  
  #### Get the probability version of the scores. ####
  probs <- exp(scores)/(1+exp(scores))
  dataset <- cbind(dataset, probs)
  colnames(dataset)[dim(dataset)[2]] <- varname2
  
  
  return(dataset)
  
}

trad_mod_match <- apply_model(model=risk_ortho_post_2010, dataset=trad_mod_match, varname1="risk_logit_ortho_post_2010", varname2="risk_prob_ortho_post_2010", era="mod")
trad_mod_match <- apply_model(model=risk_ortho_pre_2010, dataset=trad_mod_match, varname1="risk_logit_ortho_pre_2010", varname2="risk_prob_ortho_pre_2010", era="mod")
trad_mod_match <- apply_model(model=risk_gensurg_post_2010, dataset=trad_mod_match, varname1="risk_logit_gensurg_post_2010", varname2="risk_prob_gensurg_post_2010", era="mod")
trad_mod_match <- apply_model(model=risk_gensurg_pre_2010, dataset=trad_mod_match, varname1="risk_logit_gensurg_pre_2010", varname2="risk_prob_gensurg_pre_2010", era="mod")

trad_mod_match <- apply_model(model=crisk_ortho_post_2010, dataset=trad_mod_match, varname1="crisk_logit_ortho_post_2010", varname2="crisk_prob_ortho_post_2010", era="mod")
trad_mod_match <- apply_model(model=crisk_ortho_pre_2010, dataset=trad_mod_match, varname1="crisk_logit_ortho_pre_2010", varname2="crisk_prob_ortho_pre_2010", era="mod")
trad_mod_match <- apply_model(model=crisk_gensurg_post_2010, dataset=trad_mod_match, varname1="crisk_logit_gensurg_post_2010", varname2="crisk_prob_gensurg_post_2010", era="mod")
trad_mod_match <- apply_model(model=crisk_gensurg_pre_2010, dataset=trad_mod_match, varname1="crisk_logit_gensurg_pre_2010", varname2="crisk_prob_gensurg_pre_2010", era="mod")

new_exp <- apply_model(model=risk_ortho_post_2010, dataset=new_exp, varname1="risk_logit_ortho_post_2010", varname2="risk_prob_ortho_post_2010", era="mod")
new_exp <- apply_model(model=risk_ortho_pre_2010, dataset=new_exp, varname1="risk_logit_ortho_pre_2010", varname2="risk_prob_ortho_pre_2010", era="mod")
new_exp <- apply_model(model=risk_gensurg_post_2010, dataset=new_exp, varname1="risk_logit_gensurg_post_2010", varname2="risk_prob_gensurg_post_2010", era="mod")
new_exp <- apply_model(model=risk_gensurg_pre_2010, dataset=new_exp, varname1="risk_logit_gensurg_pre_2010", varname2="risk_prob_gensurg_pre_2010", era="mod")

new_exp <- apply_model(model=crisk_ortho_post_2010, dataset=new_exp, varname1="crisk_logit_ortho_post_2010", varname2="crisk_prob_ortho_post_2010", era="mod")
new_exp <- apply_model(model=crisk_ortho_pre_2010, dataset=new_exp, varname1="crisk_logit_ortho_pre_2010", varname2="crisk_prob_ortho_pre_2010", era="mod")
new_exp <- apply_model(model=crisk_gensurg_post_2010, dataset=new_exp, varname1="crisk_logit_gensurg_post_2010", varname2="crisk_prob_gensurg_post_2010", era="mod")
new_exp <- apply_model(model=crisk_gensurg_pre_2010, dataset=new_exp, varname1="crisk_logit_gensurg_pre_2010", varname2="crisk_prob_gensurg_pre_2010", era="mod")

trad_mod_match$risk_logit_mod <- ifelse(trad_mod_match$ortho_adm_mod==1 & trad_mod_match$oct_2010_plus==1, trad_mod_match$risk_logit_ortho_post_2010,
                                        ifelse(trad_mod_match$ortho_adm_mod==1 & trad_mod_match$oct_2010_plus==0, trad_mod_match$risk_logit_ortho_post_2010,
                                               ifelse(trad_mod_match$gensurg_adm_mod==1 & trad_mod_match$oct_2010_plus==1, trad_mod_match$risk_logit_gensurg_post_2010,
                                                      trad_mod_match$risk_logit_gensurg_pre_2010)))
trad_mod_match$risk_prob_mod <- ifelse(trad_mod_match$ortho_adm_mod==1 & trad_mod_match$oct_2010_plus==1, trad_mod_match$risk_prob_ortho_post_2010,
                                        ifelse(trad_mod_match$ortho_adm_mod==1 & trad_mod_match$oct_2010_plus==0, trad_mod_match$risk_prob_ortho_post_2010,
                                               ifelse(trad_mod_match$gensurg_adm_mod==1 & trad_mod_match$oct_2010_plus==1, trad_mod_match$risk_prob_gensurg_post_2010,
                                                      trad_mod_match$risk_prob_gensurg_pre_2010)))
trad_mod_match$crisk_logit_mod <- ifelse(trad_mod_match$ortho_adm_mod==1 & trad_mod_match$oct_2010_plus==1, trad_mod_match$crisk_logit_ortho_post_2010,
                                        ifelse(trad_mod_match$ortho_adm_mod==1 & trad_mod_match$oct_2010_plus==0, trad_mod_match$crisk_logit_ortho_post_2010,
                                               ifelse(trad_mod_match$gensurg_adm_mod==1 & trad_mod_match$oct_2010_plus==1, trad_mod_match$crisk_logit_gensurg_post_2010,
                                                      trad_mod_match$crisk_logit_gensurg_pre_2010)))
trad_mod_match$crisk_prob_mod <- ifelse(trad_mod_match$ortho_adm_mod==1 & trad_mod_match$oct_2010_plus==1, trad_mod_match$crisk_prob_ortho_post_2010,
                                       ifelse(trad_mod_match$ortho_adm_mod==1 & trad_mod_match$oct_2010_plus==0, trad_mod_match$crisk_prob_ortho_post_2010,
                                              ifelse(trad_mod_match$gensurg_adm_mod==1 & trad_mod_match$oct_2010_plus==1, trad_mod_match$crisk_prob_gensurg_post_2010,
                                                     trad_mod_match$crisk_prob_gensurg_pre_2010)))

new_exp$risk_logit_mod <- ifelse(new_exp$ortho_adm_mod==1 & new_exp$oct_2010_plus==1, new_exp$risk_logit_ortho_post_2010,
                                        ifelse(new_exp$ortho_adm_mod==1 & new_exp$oct_2010_plus==0, new_exp$risk_logit_ortho_post_2010,
                                               ifelse(new_exp$gensurg_adm_mod==1 & new_exp$oct_2010_plus==1, new_exp$risk_logit_gensurg_post_2010,
                                                      new_exp$risk_logit_gensurg_pre_2010)))
new_exp$risk_prob_mod <- ifelse(new_exp$ortho_adm_mod==1 & new_exp$oct_2010_plus==1, new_exp$risk_prob_ortho_post_2010,
                                       ifelse(new_exp$ortho_adm_mod==1 & new_exp$oct_2010_plus==0, new_exp$risk_prob_ortho_post_2010,
                                              ifelse(new_exp$gensurg_adm_mod==1 & new_exp$oct_2010_plus==1, new_exp$risk_prob_gensurg_post_2010,
                                                     new_exp$risk_prob_gensurg_pre_2010)))
new_exp$crisk_logit_mod <- ifelse(new_exp$ortho_adm_mod==1 & new_exp$oct_2010_plus==1, new_exp$crisk_logit_ortho_post_2010,
                                         ifelse(new_exp$ortho_adm_mod==1 & new_exp$oct_2010_plus==0, new_exp$crisk_logit_ortho_post_2010,
                                                ifelse(new_exp$gensurg_adm_mod==1 & new_exp$oct_2010_plus==1, new_exp$crisk_logit_gensurg_post_2010,
                                                       new_exp$crisk_logit_gensurg_pre_2010)))
new_exp$crisk_prob_mod <- ifelse(new_exp$ortho_adm_mod==1 & new_exp$oct_2010_plus==1, new_exp$crisk_prob_ortho_post_2010,
                                        ifelse(new_exp$ortho_adm_mod==1 & new_exp$oct_2010_plus==0, new_exp$crisk_prob_ortho_post_2010,
                                               ifelse(new_exp$gensurg_adm_mod==1 & new_exp$oct_2010_plus==1, new_exp$crisk_prob_gensurg_post_2010,
                                                      new_exp$crisk_prob_gensurg_pre_2010)))



#############################################################################
#### For each surgeon, count how many patients they saw before matching. ####





comorbidities <- c("silb_CHF","silb_PastMi","silb_PastARR","silb_ANGINA","silb_Valvulardis","silb_UnstAngina","silb_HTN","silb_DEMENT","silb_LIVER","silb_PARAPLEG","silb_RENLDYS",
                   "silb_renlfail","silb_SEIZUR","silb_STROKE","silb_AbdCancer","silb_ASTHMA","silb_CANCER","silb_CHRNLUNG","silb_DIABETES","silb_PostPulmFibr",
                   "silb_AIDS","silb_ChrPepticUlcer","silb_COAG","silb_CollagenVascular","silb_CongenitalCoag","silb_CUSHINGS","silb_GRAVES","silb_HYPOTHY","silb_SickleCell","silb_THROMB")

procs <- c("pd_access_procedure","small_bowel_resection","hernia_abdominal_open",
           "cholecystectomy", "lysis_of_adhesions", "colectomy_partial__open", "colectomy_partial__lap",
           "appendectomy", "hernia_groin_open", "large_bowel_other",  "colectomy_total__open",
           "proctectomy", "pyloroplasty", "ulcer", "mastectomy",
           "thyroidectomy_partial", "stomach_antireflux", "thyroidectomy_total", "parathyroidectomy",
           "splenectomy", "stomach_gastric_bypass_nonbariatric", "hernia_abdominal_lap", "stomach_partial_gastrectomy",
           "stomach_other", "pancreatectomy", "biliary_other", "hernia_groin_lap",
           "stomach_total_gastrectomy", "small_bowel_other",  "adrenalectomy",      "esophagectomy",
           "biliary_common_duct","liver_other", "esophagomyotomy", "bariatric",
           "liver_partial_hepatectomy","proctopexy", "thyroidectomy_substernal","colectomy_total__lap",
           "hip_repair","knee_replacement","hip_replacement","knee_replacement_revision",
           "hip_replacement_revision")

other_vars <- c("emergent_type", "admit_year","admsn_in_quarters", "comorb_sum", "age_at_admsn",
                "median_income", "pct_hsgrad", "transfer_24hrs_ip", "age_65_69", "age_70_74", "age_75_79",
                "age_80_85", "age_85_plus", "spx_major_narrow", "past_adm_6mo", "oct_2010_plus", "jan_2012_plus",
                "sex_female", "dual","binary_pct_hsgrad", "weekend")



con_vars <- c("admit_year","admsn_in_quarters","comorb_sum", "age_at_admsn",
              "median_income", "pct_hsgrad",
              "prop_prob_pat_mod_post_2010" ,"prop_log_pat_mod_post_2010",             
              "prop_prob_pat_mod_pre_2010","prop_log_pat_mod_pre_2010",
              "prop_log_pat_mod", "prop_prob_pat_mod",
              "risk_logit_mod","risk_prob_mod",
              "crisk_logit_mod","crisk_prob_mod")
bin_vars <- c(comorbidities, procs, "emergent_type", "transfer_24hrs_ip", "age_65_69", "age_70_74", "age_75_79",
              "age_80_85", "age_85_plus", "spx_major_narrow", "past_adm_6mo", "oct_2010_plus", "jan_2012_plus",
              "sex_female", "dual","binary_pct_hsgrad", "weekend", "certified_by_era")

all_vars <- c(con_vars, bin_vars)


#### Compute an overall balance table.
patient_bal_tab <- NULL
patient_bal_tab <- baltabfunc_patient(cont_vars = con_vars,
                                      binary_vars = bin_vars,
                                      before_cases=new_exp[which(new_exp$new_surgeon_flag==1),],
                                      after_cases=trad_mod_match[which(trad_mod_match$new_surgeon_flag==1),which(names(trad_mod_match) %in% names(new_exp))],
                                      before_controls = new_exp[which(new_exp$new_surgeon_flag==0),],
                                      after_controls = trad_mod_match[which(trad_mod_match$new_surgeon_flag==0),which(names(trad_mod_match) %in% names(new_exp))],
                                      bal_tab_vars = all_vars)

patient_bal_tab
write.csv(patient_bal_tab, file="/Users/sharpe/Desktop/Young Surgeons/taper_table/lv_0_overall_baltab.csv")
#### Compute era-specific balance tables. ####


mod_bal_tab <- NULL
mod_bal_tab <- baltabfunc_patient(cont_vars = con_vars,
                                  binary_vars = bin_vars,
                                  before_cases=new_exp[which(new_exp$new_surgeon_flag==1 & new_exp$era_mod==1),],
                                  after_cases=trad_mod_match[which(trad_mod_match$new_surgeon_flag==1 & trad_mod_match$era_mod==1),which(names(trad_mod_match) %in% names(new_exp))],
                                  before_controls = new_exp[which(new_exp$new_surgeon_flag==0 & new_exp$era_mod==1),],
                                  after_controls = trad_mod_match[which(trad_mod_match$new_surgeon_flag==0 & trad_mod_match$era_mod==1),which(names(trad_mod_match) %in% names(new_exp))],
                                  bal_tab_vars = all_vars)

write.csv(x=mod_bal_tab, file="/Users/sharpe/Desktop/Young Surgeons/taper_table/lv_0_mod_baltab.csv")

#### Compute surgery-specific balance tables. ####

ortho_bal_tab <- NULL
ortho_bal_tab <- baltabfunc_patient(cont_vars = con_vars,
                                    binary_vars = bin_vars,
                                    before_cases=new_exp[which(new_exp$new_surgeon_flag==1 & (new_exp$ortho_adm_trad==1 | new_exp$ortho_adm_mod==1)),],
                                    after_cases=trad_mod_match[which(trad_mod_match$new_surgeon_flag==1 & (trad_mod_match$ortho_adm_trad==1 | trad_mod_match$ortho_adm_mod==1)),which(names(trad_mod_match) %in% names(new_exp))],
                                    before_controls = new_exp[which(new_exp$new_surgeon_flag==0 & (new_exp$ortho_adm_trad==1 | new_exp$ortho_adm_mod==1)),],
                                    after_controls = trad_mod_match[which(trad_mod_match$new_surgeon_flag==0 & (trad_mod_match$ortho_adm_trad==1 | trad_mod_match$ortho_adm_mod==1)),which(names(trad_mod_match) %in% names(new_exp))],
                                    bal_tab_vars = all_vars)
write.csv(x=ortho_bal_tab, file="/Users/sharpe/Desktop/Young Surgeons/taper_table/lv_0_ortho_baltab.csv")


gensurg_bal_tab <- NULL
gensurg_bal_tab <- baltabfunc_patient(cont_vars = con_vars,
                                      binary_vars = bin_vars,
                                      before_cases=new_exp[which(new_exp$new_surgeon_flag==1 & (new_exp$gensurg_adm_trad==1 | new_exp$gensurg_adm_mod==1)),],
                                      after_cases=trad_mod_match[which(trad_mod_match$new_surgeon_flag==1 & (trad_mod_match$gensurg_adm_trad==1 | trad_mod_match$gensurg_adm_mod==1)),which(names(trad_mod_match) %in% names(new_exp))],
                                      before_controls = new_exp[which(new_exp$new_surgeon_flag==0 & (new_exp$gensurg_adm_trad==1 | new_exp$gensurg_adm_mod==1)),],
                                      after_controls = trad_mod_match[which(trad_mod_match$new_surgeon_flag==0 & (trad_mod_match$gensurg_adm_trad==1 | trad_mod_match$gensurg_adm_mod==1)),which(names(trad_mod_match) %in% names(new_exp))],
                                      bal_tab_vars = all_vars)
write.csv(x=gensurg_bal_tab, file="/Users/sharpe/Desktop/Young Surgeons/taper_table/lv_0_gensurg_baltab.csv")


